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Abstract 

The classical bond-fluctuation model (BFM) is an efficient lattice Monte Carlo algorithm for 
coarse-grained polymer chains where each monomer occupies exclusively a certain number of lat- 
tice sites. In this paper we propose a generalization of the BFM where we relax this constraint 
and allow the overlap of monomers subject to a finite energy penalty e. This is done to vary 
systematically the dimensionless compressibility g of the solution in order to investigate the influ- 
ence of density fluctuations in dense polymer melts on various static properties at constant overall 
monomer density. The compressibility is obtained directly from the low-wavevector limit of the 
static structure factor. We consider, e.g., the intrachain bond-bond correlation function, P{s), 
of two bonds separated by s monomers along the chain. It is shown that the excluded volume 
interactions are never fully screened for very long chains. If distances smaller than the thermal 
blob size are probed (s <ti g) the chains are swollen according to the classical Fixman expansion 
where, e.g., P{s) ~ g~^s~^/'^ . More importantly, the polymers behave on larger distances {s ^ g) 
like swollen chains of incompressible blobs with P{s) ~ g'^s~^/'^. 
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I. INTRODUCTION 



The bond-fluctuation model. The classical bond-fluctuation model (BFM) is an eflicient 
lattice Monte Carlo (MC) algorithm for coarse-grained polymer chains where each monomer 
occupies exclusively a certain number of lattice sites on a simple cubic lattice 
It was proposed in 1988 by Carmesin and Kremer ij] as an alternative to single-site self- 
avoiding walk models, which retains the computational efficiency of the lattice without being 
plagued by severe ergodicity problems. The key idea is to increase the size of a monomer 
which now occupies a whole unit cell of the lattice, as illustrated in Fig. [TJ The multitude 
of possible bond lengths and angles allows a better representation of the continuous-space 
behavior of real polymer solutions and melts. The BFM algorithm has been used for a 
huge range of problems addressing the generic behavior^ of long Jgolymer^ chains of very 



different molecular architectures and geometries: statics 



and dynamics |l4 . 
interfaces 



261], brushes 



self-assembly 




i6i, 0, Q 



jgels and network 



of linear chains and rings |19l 

iS 



291], thin films 



31 




10 



11 



12 



13| 



polymer blends and 



24| . glass transition [251], (co-) polymers at surfaces 
321 ] , equilibrium polymers 111 . l33l . [SJ] and general 
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38]. 



361 ]. .... For recent reviews see Refs. 
A BFM version allowing a systematic compressibility variation. As sketched in Fig. [H 
we propose here a generalization of the BFM where we relax the no-overlap constraint and 
allow the overlap of monomers subject to a finite energy penalty e. This is done to vary 
systematically the strength of density fiuctuations in dense solutions and melts to study 
their infiuence on static and dynamical properties. More specifically, we want to test the 
standard perturbation theory of weakly interacting three-dimension al p olymer melts 39] 



n n 



11 



12 



, 131 , which have been 



and to verify whether certain long-range correlations [9|, [lO, 
found recently for incompressible melts, are also present in melts with finite compressibility. 
We will see that this is indeed the case if one considers properties on scales larger than the 



screening length ^ ^ bg^^'^ of the density fiuctuations where b indicates the effective bond 




length of the chain and g the number of monomers spanning the "thermal blob" 
Interestingly, g is related to the isothermal compressibility of the solution kt [111 , 
i.e. the thermodynamic property which measures the strength of the density fluctuations. 
It can thus be determined directly experimentally or in a computer simulation from the 



2 



low-wavevector limit of the total monomer structure factor 4^, |45|, |46(] 



G{q) = Yl ■ {rn-rm))^ g = T n^p, (1) 

^mon ^^^^^^ 

with q being the wavevector, the position of monomer i, n^on = pV the total monomer 
number, p the monomer number density, V = L"^ the volume of the system and T = 1/(3 
the temperature. (Boltzmann's constant k-Q is set to unity throughout this paper.) Due 



to the definition given in Eq. ([T]), (7 is also called "dimensionless compressibility" 47j. We 
will show that the variation of the operational parameter x = e/T [48 1, characterizing the 
strength of the overlap penalty, allows us to scan g{x) over four orders of magnitude 41|. 
This puts us in the position to test various theoretical predictions which are sketched below. 

Thermodynamic properties for x ^ 1. To characterize the proposed soft BFM model 
we will first investigate thermodynamic properties such as the mean overlap energy per 
monomer e, the excess chemical potential or the (already mentioned) dimensionless 
compressibility g as functions of x. In the limit of weak overlap penalties these properties 



have been calculated long ago by Edwards theory j39| and we will compare our numerical 
results with his predictions. For later reference we postulate here the contributions to the 
free energy per monomer /(/3) relevant for this comparison 

+ i (log(p/iV) - 1) 



The first term is due to the constant intrachain self-energy discussed in Sec. IIII A[ 
due to the reference energy chosen in this study and it is not accessible experimentally 



t is 



49|. 



The second and third line of Eq. ([2]) can be readily obtained from the literature, e.g., by 



integrating the osmotic pressure given by Eq. (5.45) or Eq. (5. II. 5) of [39| with respect to 
the density p. (The first line can be considered as the integration constant with respect to 
this integration.) The second line represents the translational invariance of monodisperse 
chains of length (van't Hoff's law). Due to this contribution the compressibilities depend 
in general on as will be discussed in Sec. IIII D[ The (bare) excluded volume interaction 
between the monomers is accounted for by the first term in the third line with v{x) being 
the second virial coefficient of a solution of unconnected monomers. The underlined term 



3 



represents the leading correction to the previous term due to the fact that the monomers 
are connected by bonds summing over the density fluctuations to quadratic order. As one 
expects {4^, the corresponding correlations of the density fluctuations reduce the free eriergy 
by about one k^T per thermal blob of volume Consistently with Refs. 
correlation length ^ of the density fluctuations has been defined here as 
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39, 



43j the 



^^{x) = b^{x)g{x)/12 ^ b'^{x)/{l2v{x)p). 



(3) 



Note that g{x) and b{x) refer respectively to the dimensionless compressibility and the effec- 



41 



50| . The density fluctuation contribution 



five bond length of asymptotically long chains 
to the free energy will be demonstrated numerically from the scaling of the specific heat cy 
with respect to overlap strength x and density p (Sec. 1111 B|) . We note finally that Eqs. ([2]) 
and ([3]) are supposed to apply as long as 



G,{x) = 



1 



'g{x) h'^{x)p 



«1, 



(4) 



with the Ginzburg parameter being the small parameter of the perturbation theory 51|. 
This restricts — as we shall see — the related predictions to overlap penalties with x ^ 1. 

Conformational properties of asymptotically long chains for all x. This paper aims ulti- 
mately to characterize an important intrachain conformational property, the size of chain seg- 



ments of arc-length s in very long chains where finite- effects may be neglected 
ically, we will investigate the mean-square end-to-end distance -R^(s) = ((r. 



52|. Specif- 



m=n+s 



where the average is performed over all possible pairs of monomers (n, m = n + s). The 
total chain end-to-end distance is Rc{N) = R{s = N — 1). If appropriately plotted 12| . 



the segment size R(s) will allow us to obtain an estimation of the effective bond length b{x) 
by extrapolation 4l|- Our numerical results will be compared with an analytical predic- 
tion obtained by a standard one-loop perturbation calculation very similar to the analytic 



lol. 
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Details of the calculation for gen- 
53|. Focusing in this paper mainly on 



results already presented elsewhere 
eral soft melts will be presented in a future paper 
computational issues, we only quote here the "key relation" put to the test 



1 - 



1 



9 



1/2 



u 



u 



(l-e'"erfc(y2^))j 



(5) 



where we have used u = s/g{x) for the reduced curvilinear arc-length and erfc(x) for the 
complementary error function ^]. The correction to Gaussianity expressed by the r.h.s. 



terms of the key prediction is always positive and corresponds to a weak swelling of the 
chain. The "swelhng coefficient" Cg = \p2AfTx^ j pi? 12| measures the strength of the effect. 
While Eq. ((21) fails for large x where the Ginzburg parameter Gz becomes of order unity, 
Eq. ([5]) is supposed to hold for all x provided that the considered segment length s is large 
enough. As explained in Sec. II of Ref. pL2], the relevant small parameter of the perturbation 
theory is here the so-called "correlation hole potential" of the chains 



u 



» ~ G^^/gfs ^ l/^/s < 1 for s > 1 



(6) 



rather than G^. For a discussion of chain end effects in incompressible melts see Ref. 121]. 

Check of limiting behavior. Although no formal derivation of Eq. (jS]) is given, we verify 
briefly whether the suggested key prediction is reasonable by discussing the limits of large 
and small reduced arc-length. We remind first that according to Flory's "ideality hypothesis" 



55| . polymer chains in the melt are thought to display Gaussian statistics for segment 



sizes somewhat larger than the persistence length |39l . 1401 . |55| which implies that the r.h.s. 
of Eq. (I5l) is traditionally assumed to vanish (exponentially) for finite s. At variance to this, 
the expansion of the complementary error function 5J] for u = s/g ^ 1 yields 

1 _ ^ ~ ^ (7) 

suggesting in fact that corrections to Gaussianity must be taken into account for all finite 
s. As one expects, however, Gaussianity is still recovered in three dimensions if s — oo. 
(Incidentally, this does not hold for effectively two-dimensional melts as may be seen from 



the discussion of ultrathin polymer films in Refs. |32l. |42|. |56| . ) Most remarkably, the explicit 
compressibility dependence drops out for large u for all g (x) provided that the chains are long 
enough, such that g <^ s <^ N. Eq. ([7]) is precisely the relation which has been discussed in 



detail both theoretically and numerically 
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12j for highly incompressible melts with 



full excluded volume interactions (x = oo). Hence, this is the expected limiting behavior if 
the polymer chains are renormalized in terms of an incompressible packing of thermal blobs 
with g monomers. In the opposite limit of small reduced arc- lengths, u = s/g <^ 1, the 
expansion of Eq. ([5]) yields 

. R'^(s) f R 4 



V~9 



(8) 



This is consistent with the classical expansion result of the chain size in terms of the "Fixman 
parameter" z ~ Vy/s/b^ ~ V^/ {y/gb^p) [391, |47| since the deviations from ideality expressed 



by the last term of Eq. ([8]) become then proportional to —Cs\Ju/g ~ — z, in agreement with 
the leading correction term for the total chain size Re{N) presented in textbooks 57|. 



Outline. In Section [TTl the algorithm is introduced and some technical details are dis- 
cussed. First we summarize the classical BFM without monomer overlap (Sec. Ill AD . which 
has been used in previous work providing the start configuration for the present study, 
and introduce then its generalization with finite overlap penalty (Sec. IIIBI) . The cen- 
tral Section IIIII presents our numerical results starting with the thermodynamic properties 
(Sees. UTTAIIIII DI) ■ in particular the dimensionless compressibility g{x). We characterize then 
(Sees, nil EtHTTGl) various intrachain properties, as for instance the effective bond length b{x) 
or the bond-bond correlation function P{s), as functions of g{x). A synopsis of our results 
is presented in Section [IVl 



II. ALGORITHM AND TECHNICAL DETAILS 

A. Bond fluctuation model without monomer overlap 

A lattice Monte Carlo scheme for topology conserving polymers. We have used the three 



di.e.io„a, ve.io„ of Bond fluc.,a.io„ .ode, ^ S fl M, ^ whe. each effective 

coarse-grained monomer is represented by a cube of eight adjacent sites on a simple cubic 
lattice, as illustrated in Fig. [H Even the partial overlap of monomers is forbidden. The 
lattice constant a is naturally chosen as the unit length. Polymers of length consist of 
cubes connected by — 1 bonds, as shown in the sketch for A^ = 3. These bonds are taken 
from the set 

{P(2, 0, 0), P(2, 1, 0), P(2, 1, 1), P(2, 2, 1), P(3, 0, 0), P(3, 1, 0)} (9) 

of allowed bond vectors, where P stands for all permutations and sign combinations of 
coordinates. This corresponds to 108 different bond vectors of 5 possible bond lengths (2, 
\/5, ^/6, 3, vTO) and 100 possible angles between consecutive bonds. The smallest 13 angles 
do not appear for the classical BFM without monomer overlap, because excluded volume 
forbids sharp backfolding of bonds. If only local Monte Carlo (MC) moves of the monomers 
to the six nearest neighbor sites are performed — called "L06" moves [l2| — this set of 
vectors ensures automatically that polymer chains cannot cross. Topological constraints. 



e.g. in ring polymers 



19| , hence are conserved and the polymer dynamics may be expected 



to be of reptation type {4, 12, Q, Q. It is this fact which has originally motivated the 
choice of allowed bonds. We keep it for consistency with previous work although the non- 
crossing constraint is irrelevant for the present study. Note that the classical BFM algorithm 
with L06 moves is strictly speaking not ergodic, since some (thermodynamically irrelevant) 
configurations may be easily constructed which are not accessible starting from an initial 
configuration of stretched linear chains 37l |. 

Obtaining athermal start configurations with topology violating moves. The algorithm is 
up to now athermal and the only control parameters are the monomer density p and the 
chain size N . Melt conditions are realized for p = 0.5/8, where half of the lattice sites are 
occupied {4]. We use (if not stated otherwise) periodic simulation boxes of linear dimension 
L = 256 which contain n^on = pL^ = 2^'^ ~ 10^ monomers. This large system size allows to 
eliminate finite-size effects even for the longest chain lengths studied. Our simulations have 



61 



62 



been carried out by a mixture of local, slithering snake [17|, |59, 6^, and double-bridging MC 



moves [37|, l6l|, [6^, l63|] which allow us to equilibrate polymer melts with chain lengths up 
to = 8192. Instead of the more realistic but very slow L06 dynamical scheme we make 
so-called "L26" jump attempts to the 26 sites of the cube surrounding the current monomer 
position. This permits the crossing of chains which dramatically speeds up the dynamics, 
expecially for long chains {N > 512). Details of the equilibration procedure and possible 



caveats are discussed in Ref. 



12|. We stress that if these topology violating MC moves are 



included all possible configurations become accessible, i.e. the BFM becomes fully ergodic. 
Concerning the static properties ergodic and non-ergodic BFM versions are, however, prac- 
tically equivalent. This has been confirmed by comparing various static properties and by 
counting the number of monomers which become "blocked" once one returns to the original 
local L06 moves. 



B. Bond fluctuation model with finite excluded volume penalty 

The model Hamiltonian. Fig. [1] shows how finite energy penalties are introduced. The 
overlap of two cube corners on one lattice site (A^ov = 1) corresponds to an energy cost 
of e/8, the full overlap of two monomers {No-^ = 8) to an energy e. More generally, with 
A^ov being the total number of interacting cube corners the total interaction energy of a 
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configuration is 



E 



(10) 



With the energies of the final (Ef) and the initial configurations [E]) we accept the MC move 
according to the Metropolis criterion 37|, l6J] with probability min(l, exp[— (E'f — Ei)/T]). 
We set arbitrarily e = 1 and vary the ratio x = e/T 48|| starting from x = oo corresponding 
to the athermal classical BFM and systematically increase the temperature T as shown in 
Table [I for unconnected monomers (A^ = 1). 

Second virial coefficient. To illustrate this interaction we indicate already here the sec- 
ond virial of an imperfect gas of unconnected monomers, v = J d6{l — e~^^^'^^'^), which is 
shown below to be useful for roughly characterizing the effective strength of the potential. S 
stands for a possible lattice vector between the centers of two interacting cubes. It is easy to 
see that there are 8 vectors corresponding to A'ov = 1 (as in Fig. [1]), 12 to N^-^ = 2 (overlap 
of two cube corners), 6 to N^^ = 4 (overlap of two faces), and 1 to No-v = 8 (full overlap). 
This leads to a second virial 



V 



8x (l-exp(-a;/8)) + 12x (l-exp(-a;/4))+6x (l-exp(-a;/2)) + lx (l-exp(-x)) (11) 



given in units of the lattice cube volume a^. We note that the second virial becomes constant, 
V = 27, in the low temperature limit (x ^ 1) corresponding to the classical athermal BFM 
result 



2[. In the opposite high temperature limit {x <C 1) it decays as 

27 o 



iX 



16^ 



-X 



(12) 



Implementation of the algorithm. We briefly explain the implementation of BFM chains 



with the soft overlap Hamiltonian, Eq. ( ITOi) . Following [J, |l4, Il6l . Il7l | it is convenient 
to keep one list of the monomer positions in absolute space (since we are also interested 
in dynamical properties) and one for the corresponding lattice positions. We identify each 
of the 108 allowed bonds of the set [Eq. (Q] with a unique bond index and keep a list of 
these indices. This allows to verify rapidly whether an attempted bond vector is acceptable. 
Since the bond index (being less than 128) can be encoded as a byte, this compressed list is 
stored when we write down the configuration ensemble for further analysis. Additional lists 
allow to handle efficiently the periodic boundary conditions, the change of a bond index for 
a given monomer move and the local interactions relevant for local L06 or L26 moves. 
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Since the soft overlap Hamiltonian allows the occupation of a lattice site by more than 
one monomer, it is not possible to use a compact boolean occupation lattice (corresponding 
to a spin cr = and cr = 1 for an empty or filled lattice site) or an integer lattice filled with 
the monomer indices as in previous implementations. Instead we have mapped Eq. flTU]) onto 
a Potts spin model [6^] 

r <5 

with constant monomer number TT-mon = y^.^ cr(^) = L^P using a Wigner-Seitz representation 
of the cubic lattice following Miiller 2l|, |38|] . In this representation an integer spin variable 
o"(r) counts the number of BFM monomers (cr = 0,1,2,.. .) with cubes centered at a Wigner- 
Seitz lattice position r. In other words, each cube is not represented by 8 lattice entries for 
the cube corners, but just by one for its center. Since we have now to compute the interaction 
between cube centers instead of cube corners, the coupling constant J characterizing the 
interaction between two spins depends only on the relative distance 6: 

f 

1/8 if 5 = P(l, 1, 1) (cube corners), 
1/4 if (5 = P(l, 1, 0) (cube edges), 
J((5) =e{ 1/2 if (5 = P(l, 0, 0) (cube faces), (14) 
1 if 5 = P(0, 0, 0) (full overlap), 
otherwise. 

Since the interaction is still short-ranged and the values of J are readily tabulated, this 
remains an efficient rendering of the monomer interactions. Note that the first term on the 
r.h.s. of Eq. (fT3l) contains a constant self- interaction contribution of the rimon monomers 
with themselves for (5 = 0, which is substracted by the second term j65 |. 

Equilibration and system properties of high-molecular melts. As already stated we have 
used as start configurations the equilibrated BFM configurations without monomer overlap 
from our computationally much more expensive previous studies obtained with topology 
violating moves pL2i]. Since the soft BFM simulations are also ergodic, these are the relevant 
reference configurations. Starting with these configurations we increase the temperature. 
As one may expect, the configurational properties essentially are found unchanged for low 
temperatures (x ^ 5). Local L26 moves need to be added to global slithering-snakes moves 
for a; > 1. Otherwise the pure snake motion will become ineffective as it is well known from 



a previous study of the snake dynamics without overlap 17|] . Simple slithering-snake moves 
without additional local moves are sufficient for smaller x. We have crosschecked our results 
in this regime for N = 2048 and = 8192 using boxes of linear size L = 512 by starting our 
simulations with Gaussian chains at a; = and increasing then the penalty. Tables [ITl and IIIII 
present some system properties obtained for our reference density p = 0.5/8. Averages are 
performed over all chains and at least 100 configurations. Table HTl summarizes the properties 
extrapolated for asymptotically long chains. Similar information is given in Table IIIII for a 
constant penalty function of chain length A^. Density effects have been studied 

only briefly for chains of length A^ = 8192 and weak overlap penalties (x ^ 1). This has 
been done to investigate the intrachain contributions to the mean energy. We begin the 
discussion of our numerical results by addressing this issue. 



III. COMPUTATIONAL RESULTS 



A. The mean overlap energy 



Qualitative behavior. From the numerical point of view the simplest thermodynamic 
property to be investigated here is the mean interaction energy per monomer, e = (E) /n^^on, 
due to the Hamiltonian, Eq. (ITU]) . Fig. [2] presents the dimensionless energy y = e/e as a 



function of the reduced overlap penalty x = e/T = e(5 48|] for melts at monomer number 



density p = 0.5/8 and various chain lengths A^ as indicated. We increase the temperature T 



from the right to the left starting with configurations obtained recently [l2j for the classical 
athermal BFM. As one expects, the interaction energy increases exponentially for small T 
and levels off for large T where chains and their monomers freely overlap (x ^ 1). The data 
for unconnected beads (A^ = 1) represented by the ffiled spheres and polymer chains (A^ ^ 1) 
are broadly speaking similar, especially for large overlap penalties, x > 1. Interestingly, the 
mean energy of polymer melts increases more strongly for x <^ 1 as can be seen better from 
the log-linear data representation chosen in the inset of Fig. [2l Also shown in the inset is 
the mean intrachain self-energy per chain monomer Cgeif (filled triangles) obtained for the 
largest chain length A^ available for a given x. As can be seen also from Table [Tll about 
half of the energy of polymer melts for all x is due to these intrachain interactions. For 
X -C 1 the self-energy becomes eseii/s ~ 0.18 which is exactly the observed energy difference 



10 



between polymer and bead systems. 

Second virial contribution. Before addressing this point let us first consider the energy 
of unconnected soft BFM beads for which the second virial coefficient v{x) has been give n 



in Eq. flTT]) . Since e = dp{(3f{(3)) the mean energy per bead becomes to leading order 44l | 



1 dv{oc\ p 

y ^ -p — - — = - (exp(— x/8) + 3 exp(— x/4) + 3 exp(— x/2) + exp(— x)) (15) 

2 ox 2 

corresponding to the first term in the third line of Eq. ([2]). Eq. f|T5l) is represented by the 
dashed line in Fig. [2l It corresponds to a Arrhenius behavior with y ^ pexp(— x/8)/2 
(dash-dotted line) in the low temperature region and, as expected, to y ^ |8p = 4p for 
large temperatures. This simple formula predicts well the bead data over the entire range 
of X (underestimating slightly the mean energy at x ~ 10) and yields also a remarkable fit 
for polymer chains with larger overlap penalties. 

Self- energy in the high temperature limit. The above-mentioned energy difference be- 
tween polymer chains and beads for x -C 1 has been accounted for by the first free energy 
contribution indicated in Eq. ([2]). This contribution is further investigated in Fig. [3] pre- 
senting data for such a high temperature (x = 0.001) that the entropy dominates essentially 
all conformational properties. The self-energy of a chain is thus given by the probability 
p{s, S) that a random walk of s BFM bonds returns to a relative position 6 with respect to 
a reference monomer at r. The self-energy per monomer is then 

<5 s=2 

where the first sum runs over all positions with non-vanishing coupling constant J{S) as 
defined in Eq. ( Till) . The probability p{s, 6) and the weights J{6)p{s, 6) can be tabulated in 
principle for small s. Since the return probability decreases strongly with s, these model- 
specific small-s values dominate the integral, Eq. f|T6|) . As can be seen from the inset of 
Fig. [3] for single chains (corresponding to an overall density p = 0) the self-energy per 
monomer becomes Cseif ~ 0.18£: for large A^. The weak chain length dependence visible in 
the panel stems from the upper integration boundary over the Gaussian return probability 
which leads to a chain length correction linear in t = — 1. This is indicated by the 

bold line presented in the panel. Also shown in the panel are energies for our reference 
density p = 0.5/8. They are shifted vertically by the mean field energy 4p assuming that 
densities fiuctuations of different chains do not couple. The main panel presents the mean 
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energy e and the mean self-energy Cseif as functions of the density p for chains of length 
= 8192. The self-energy (triangles) stays essentially density-independent. The total 

interaction energy sums over the self-energy and mean-field energy contributions as shown 

by the dashed line. The self-energy contribution can only be neglected for very large densities 

corresponding to volume fractions larger than unity. 

Temperature dependence in the high temperature limit. Summarizing Eqs. ([2]), ([3]) and 

f|T2|) the energy per bead should scale to leading order in x as 



where the two x-independent contributions have already been discussed above. The under- 
lined term stems from the density fluctuation contribution for long polymer chains predicted 
by Edwards [39] indicated in Eq. ([2]). Here we have approximated the effective bond length 
h{x) by the mean-squared bond length / ~ . As can be seen from Table [ITl this approxi- 
mation (further discussed in Sec. IIIIF|) is justified for x <^ 1. Eq. f|T7j) is indicated by the 
bold line in the inset of Fig. [2J It yields a reasonable description of the temperature depen- 
dence of the mean energy for small x. Since the energy is dominated by the two constant 
contributions to Eq. ( JT71) for x < 0.001 and since higher expansion terms become relevant 
for X > 0.1, the predicted v^-decay corresponds unfortunately only to a small x— regime. 
In order to show that it is indeed the density fluctuation term which dominates the temper- 
ature dependence for x -C 1 we will consider in the next paragraph the specific heat cy, i.e. 
the second derivative of the free energy with respect to jS. 

B. Energy fluctuations 

Specific heat. The fluctuations of the interaction energy are addressed in Fig. H] display- 
ing the enthalpic contribution to the specific heat per monomer, cy = — /3^5^(/3/(/5)) = 
^ {{E"^) — (E)'^) /^mon [44]. Using again the second virial of soft BFM beads, Eq. (fTTj) . one 
obtains the simple estimate for the specific heat 



represented by the dashed line. In the large-x limit, this corresponds to the exponential 
decay, cy ~ px^ exp(— x/8)/16, indicated by the dash-dotted line. For barely interacting 




(17) 
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beads (x <^ 1), Eq. (fT8|) yields a power-law limiting behavior, cy ~ ffp^;^ ~ As can 

be seen from the plot, Eq. f|T8|) predicts the energy fluctuations of BFM beads for essentially 
all X, slightly underestimating again the maximum of cy at a; ~ 10. Since the speciflc heat 
for beads and polymer chains is similar for x > 1, the virial formula is also good for polymer 
chains in this limit. 

High temperature limit for polymer melts. Strong chain length effects are, however, visi- 
ble for high temperatures (x ^ 1) where the speciflc heat is found to increase monotonously 
with A^. This can better be seen from the inset where the speciflc heat is plotted as a func- 
tion of the reduced chain length u = N/g with g being the dimensionless compressibility 
determined below in Sec. IIII D[ (Since e and cy correspond to different derivatives of the 
free energy / with respect to f3, there is obviously no thermodynamic inconsistency in the 
flnding that cy reveals much larger chain length effects than e.) For large chains with u ^ 1 
this increase levels off at a chain length independent envelope 

as anticipated by the density fluctuation contribution predicted in Eq. ([21). In contrast 
to Eq. f|T7j) for the mean energy the density fluctuation term does now correspond to the 
leading contribution to the numerically measured property. This increases the range where 
the density fluctuation contribution can be demonstrated to over three decades in x. Eq. f[T9|) 
is indicated by the bold lines in the main panel and the inset of Fig. |H 

Scaling with chain length N. We have still to clarify the scaling observed for u = 
N/g{x) ^ 1 in the inset of Fig. HI Chains which are smaller than the thermal blob {u <^ 1) 
behave obviously as random walks and the density fluctuations decouple from the interaction 
strength x. Due to this factorization the speciflc heat for these short chains must scale as 
x^, just as for beads. This is shown in the main flgure for = 16 (thin solid line). Consis- 
tency with Eq. f|T9l) implies the scaling cy(x, A^) ~ x^/^p^/^/i(u) with h{u) being a universal 
function scaling as h{u) ~ in the large-M limit. Since cy ~ x^ and g{x) ~ l/{xp) (as 
shown in Sec. IIIIDI) for m <^ 1, it follows that h(u) ~ m^^^ as conflrmed by the dashed line 
indicated in the inset. Hence, cy ~ px'^N^/'^ for n ^ 1. 
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C. Chemical potential 



Scaling of the chemical potential. Fig. [5] presents the excess chemical potential per 
monomer, y = fi^x/TN, obtained using thermodynamic integration (as explained below) 
for three chain lengths = 1, 64, and 2048 as functions of the overlap penalty x = e/T. As 
one expects, the chemical potential increases first linearly with x and then levels off. Chain 



length effects are again small on the logarithmic scale chosen in the plot [66|. For large x the 
chemical potential becomes slightly larger for beads {y ~ 2.64) than for long chains where 
y ^ 2.1 (dash-dotted line). That the chemical potential of polymer chains is reduced com- 
pared to melts of unconnected beads is of course expected for all x due to the (effectively) 
attractive bond potential. For weak interactions this reduction should be described by the 
density fluctuation contribution to the free energy [Eq. ([2])] which corresponds to an excess 
chemical potential 

y = ^^^v{x)p^l- -^^^ + •••j forx«l (20) 

with v{x) being the second virial of unconnected beads. The dashed line in Fig.Opresents the 
leading contribution v{x)p for unconnected beads, the bold line in addition the underlined 
connectivity contribution given in Eq. (l20l) . It turns out that the simple second virial 
approximation provides a much better fit of the data over the entire x-range than the 
full prediction. (The weak underestimation of the chemical potential for x > 10 must be 
attributed to higher order correlations relevant in this limit.) That the density fluctuation 
contribution overestimates the reduction of the chemical potential for x > 1 is in agreement 
with Eq. (jl]) and the Ginzburg parameters indicated in Table HIl Hence, Eq. (1201) in principle 
can be tested only for x <^ 1. Unfortunately, in this limit the relative correction, scaling 



as Y^x/p, becomes too small to allow a fair test of the theory. A systematic variation of 
the density and an improved numerical accuracy of the chemical potentials measured are 
warranted to achieve this goal. 

Thermodynamic integration. We now explain how the data of Fig. [5] have been obtained 
numerically. The simple insertion method due to Widom j^] obviously becomes rapidly 
inefficient with increasing x, especially for longer chains. Slightly generalizing the method 
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suggested in 5|, |6| we therefore have performed a thermodynamic integration 45 1 



/?/iex = (21) 
Jxie) ^ 

over discrete values of the interaction affinity A = exp(— esg/5/8) characterizing the excluded 
volume interaction of a ghost (g) chain that is gradually inserted into an equilibrated system 
(s). (Nsg) refers to the mean number of lattice sites where system and ghost monomer cube 
corners overlap at a given interaction A. Generalizing the Potts spin mapping, Eq. (|T3l) . of 
the excluded volume interactions for homopolymers presented above, we use now two spin 
lattices, ffsir) describing (as before) the interaction of the system monomers and o"g(r") the 
ghost chain. The spin lattices are kept at the same temperature and are both characterized 
by the same (arbitrary) energy e = 1 which has to be paid for a complete overlap of two 
system monomers or two ghost monomers. The interaction of both spins is described by 

AE,, = ^s(r) '^sg(^)^g(^ + S) (22) 

with coupling constants Jsg(^) ~ £sg defined as in Eq. (HM taken apart the energy parameter 
e which is replaced by the tuneable interaction energy Ssg- Starting with decoupled system 
and ghost configurations at ^sg = 0, i.e. A = 1, we gradually increase the interaction 
parameter up to Egg = e, i.e. \{e) = exp(— e/J/8), always keeping the coupled system at 
equilibrium. Monitoring the distribution of the number Nsg of overlaps between system and 
ghost cube corners we use multihistogram methods as described in js, to improve the 
precision of the integral. 

The mean overlap number (A^^sg) (devided by 8A^) is shown in the inset of Fig. [5] as a 
function of A for = 2048 and two inverse temperatures x = 3 and x = 100. Starting from 
A = 1 the overlap number decreases monotonously with increasing coupling between system 
and ghost monomers, i.e. decreasing A. Interestingly, a power law behavior {Nsg) /N ^ A^^^ 
is found empirically for large x ^ 10 (dashed line). Fitting this power law and integrating 
then analytically Eq. (1211) provides a useful crosscheck of the numerical integration using the 
multihistogram analysis. This is a technically important finding, since the multihistogram 
analysis requires overlapping distributions of A^gg and hence much more equilibrated inter- 
mediate values A as indicated for x = 100. A detailed explanation for the observed power 
law still is missing, but it is presumably due to the systematic screening of the long range 
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correlations of the ghost chain which is swollen at A = 1 becoming more and more Gaussian 
as it feels the compression due to the surrounding host chains 4G|. 



D. The compressibility 

Compressibility g{x,N) and excess compressibility g^xix). To test the key relation 
Eq. announced in the Introduction we need accurate values for the dimensionless com- 
pressibilities g{x) = limAT^oo g{x, N) of asymptotically long chains. As suggested by Eq. ([T]), 
we compute first the dimensionless compressibility g{x,N) = limq^QG{q) from the low- 
q limit of the total monomer structure factor for different overlap penalties x and chain 
lengths (see below for details). These raw data are presented in Fig. [6] as a function of x. 
As one expects, g{x, N) decreases monotonously with overlap strength x. In contrast to the 
thermodynamic integration performed for the chemical potential [Eq. 0211) ]. the structure 
factor measures the complete compressibility, not just the excess contribution. The strong 
A^-dependence visible in the plot thus is expected from the translational entropy of the 
chains. As can be seen, e.g., from Eq. ([2]) or from the virial expansion of polymer solutions 



40], the compressibility can be written in general as 



1 _ ^ d\(3f{(3)p) _ 1 1 



g{x,N) ^ dp^ N g,,{x,N) 

for all X with gex{x, N) being the^ excess contribution to the compressibility which may, at 
least in principle, depend on A^ [66]. As can be seen from the inset of Fig. [6l all compressibil- 
ities collapse, however, on one A^-independent master curve if one plots l/g{x, N) — 1/N as 
a function of x, even the compressibilities obtained for unconnected beads (A^ = 1). Within 
numerical accuracy the A^-dependence observed for g{x,N) can therefore be attributed to 
the trivial osmotic contribution and the excess compressibility gcx ~ is thus identical 
to the compressibility g{x) of asymptotically long chains. The bold line indicated in the 
inset presents the best values of g{x) summarized in Table HTl These values have been ob- 
tained from the excess compressibilities for the largest chain length available for x > 0.001. 
A precise numerical determination of gcx{x) becomes impossible for even smaller overlap 
penalties. We thus have used for the smallest x-values sampled the theoretical prediction 

-Wp|1- forx«l (24) 



g{x) \ 2n b{xYp 
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due to the free energy [Eq. ([2])] postulated in the Introduction. The prefactor v{x)p rep- 
resenting the bare monomer interaction is indicated by the dashed hne in the main panel 
of Fig. O Hence, g{x) ~ l/(8xp) for weak interactions, i.e. the compressibility increases 
linearly with temperature. The underlined term is the leading correction due to the density 
fluctuation contribution to the free energy. It implies that the excess compressibilities for 
polymer melts and unconnected beads cannot be completely identical. However, as before 
for the chemical potential, the difference is far too small to be measurable in the limit where 
Eq. (!24|) applies. Although this result is unfortunate from the theoretical point of view, the 
data collapse observed in the inset suggests that it is acceptable to numerically estimate the 
long chain compressibility g{x) by computing the structure factors of rather short chains. 

Total monomer structure factor. We now turn to the total structure factor G{q) shown in 
Fig. [7] to explain how the compressibilites g{x, N) have been obtained. Only chains of length 

= 2048 are presented for clarity. The total monomer structure factor is obtained by com- 
puting G{q) = ^^-^ ([Xln '^'^^(q ■ T„)]^ + [Xln ^^^(^ ■ '^nYf') wherc the sums run over all the 
nmon monomers of the box and the wavevectors are commensurate with the cubic box of linear 
length L. Since the smallest possible wavevector is 27r/L, it thus is important to have large 
box sizes to scan over a sufficiently important g-range allowing a reasonable determination 
of (^(x, A^). Note that around and above g ~ 2 monomer structure and lattice effects become 
important. Since only smaller wavevectors are of interest if one is interested in universal 
physical behavior, we will focus below on wavevectors g < 1. For comparison, we have also in- 
cluded the single chain form factor F{q) = ([J2n=i cos(qr ■ + [J2n=i ■ T'n)]^^ 39| 
for X = 0.001 (bold line). Note that the qualitative shape of F{q) — decaying monotonously 
with q from its maximum value F{q = 0) = N — depends very little on the temperature 



(not shown). We remind that the "random phase approximation" (RPA) formula 39l. |40| 



^ +7^ (25) 



G{q) F{q) g{x) 

relates the total structure factor to the measured form factor. Eq. fl25|) is of course consistent 
with Eq. (!23|) in the q ^ limit. It allows to directly fit for the excess compressibility 
gex{x) using the measured structure factor G{q) and form factor F{q), at least in the ar- 
range where the RPA approximation applies. As may be seen from the figure, G{q) indeed 
decreases systematically with x, i.e. with decreasing g{x). For large temperatures (x < 3) 
it also decays monotonously with q, again in agreement with Eq. (|25|) . Interestingly, this 
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becomes qualitatively different for larger excluded volume interactions {x > 3) where the 
total structure factor is essentially constant (in double-logarithmic coordinates), very weakly 
increasing monotonously with q. The RPA formula apparently does not apply in this limit 
in agreement with Eq. (jl]). Fortunately, this is of no concern for our main purpose — to 
compute g{x) — since in precisely this limit the compressibility is readily obtained from 
a broad plateau (even for much smaller boxes) which in addition becomes chain length 
independent, as we have already seen from the inset of Fiq. iQ. Using boxes with L = 256 
it is possible to directly measure the plateau values for x < 0.3. For smaller x we have 
simulated boxes with L = 512 containing nmon ~ 8.4 ■ 10^ monomers and corresponding 
to a smallest wavevector q ^ 0.01. This box size becomes again insufficient for the largest 
temperatures we have simulated, as shown in Fig. [7] for x = 0.001 (dashed line). It is for 
these values where the RPA formula, Eq. fl25l) . allowing to fit the deviation from the (barely 
visible) plateau, has been particulary useful. 

Approximated RPA formula. We finally note that in the intermediate wavevector regime 
(where q corresponds to distances much smaller than the radius of gyration and much larger 
than the monomer size) the general RPA Eq. ( 125|) may be written as 

^ ^ + ^bHx)q' = 4- {I + iqO') (26) 



Giq) g{x) 12 ^ g{x) 

which justifies the definition given in Eq. ([3]) for the correlation length of t 
tions ^. Here we have used that the form factor becomes F[q) k, 12/6^q^ 



le density fluctua- 



391 ]. This assumes 



that corrections to Gaussian chain statistics may be ignored |10l . Ill| and that finite chain 
size effects are negligible. From the numerical point of view the approximated RPA Eq. fl26|) 
has the disadvantage that the effective bond length h{x) needs to be determined first. As 
shown in Fig. [HI it has the advantage that it allows for an additional test of the values of 
g{x) and h{x) indicated in Table HTl The main panel presents the rescaled structure factor 
G{q)/g{x) for chains of length = 8192 as a function oi Q = q^ with ^ being obtained 
from g{x) using Eq. ([3]). All data collapse on the master curve 1/(1 + Q^) indicated by 
the bold line provided that the wavevector q remains sufficiently small and no local physics 
is probed. That the used compressibilities are accurate is emphasized further in the inset 
where g{x)/G{q) — 1 is plotted as a function of using only sufficiently small wavevectors 
q. According to Eq. ( l26l) all data should collapse on the bisection line in double-logarithmic 
coordinates if the correct compressibilities are used. This is indeed the case. Unfortunately, 
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even this rather precise method still does not allow to demonstrate the density fluctuation 
contribution in Eq. fl24l) for a; ^ 1 since the same scaling collapse is obtained for the simple 
choice 1/ g{x) = v{x)p. Please note the weak deviations visible for x = 1 which is due to the 
breakdown of the RPA formula for large x mentioned above. 



E. Bond properties 

Up to now we focused on some thermodynamic features of the soft BFM model, i.e. on 
large-scale properties. Turning to configurational issues we begin by characterizing local- 
scale features of the algorithm. (Readers only interested in universal properties may wish 
to skip this paragraph.) 

Mean bond length. By definition of our version of the BFM algorithm the bond length 
is allowed to fluctuate strongly between 2 and vTo. One expects that switching on the 
excluded volume interaction, i.e. decreasing the temperature, will suppress large bonds due 
to the increasing pressure. The mean bond length commonly is characterized by the root- 
mean-square length / = (Z^)^^^. (Other moments would yield similar results.) The mean 
bond length rapidly becomes (A^ > 20) chain length independent ^. As a function of 
overlap penalty x the bond length shows a monotonous decay between x ~ 3 and x ~ 20 as 
can be seen from Fig. [9l As other local properties, / becomes constant in the small-x and 
large-x limits (dashed horizontal lines). The value l{x = 0) = 2.718 gives the lower bound 
for the effective bond length b{x) of asymptotically long chains (stars) obtained below. 

Mean bond angle and local chain rigidity. Defining the bond angle 6 between two sub- 
sequent bonds by the scalar product cos(^^) = e„ • e„+i of the normalized bond vectors 
Si = li/\li\, the local chain rigidity may be characterized by (6) and (cos(0)). Note that (6) 
and {cos{9)) can be regarded as chain length independent, just as the mean bond length. 
As can be seen from Table [Tll the local rigidity is negligible for x <^ 1, i.e. (6) ~ 90° and 
(cos(^^))) ~ due to the symmetry of the distribution p{6) with respect to 90°. The rigidity 
then increases around x ~ 1 and becomes constant again for large x where (6) ^ 82.2° and 
{cos{6)) ^ 0.106. The increase of the local rigidity for larger excluded volume interactions is 
of course expected due to the suppression of immediate backfoldings corresponding to bond 



angles 6 > 143° 15|. The distribution p{6) therefore becomes lopsided towards smaller 6 



(not shown). It is well known 39| that for chains characterized by the "freely rotating" (FR) 
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chain model such a local rigidity would lead to an effective bond length b{x) = /(x)^/cfr 
with cfr = (1 + (cos(^^)))/(l — (cos(6'))). This simple model, indicated by the crosses in 
Fig. [HI yields a qualitatively reasonable trend (monotonous increase of the effective bond 
length at x ^ 1) but fails to fit the directly measured effective bond lengths quantitatively. 



F. Chain and segment size 

Total chain size Re{N). One way to characterize the total chain size is to measure 
the second moment of the chain end-to-end distance Rl{x,N) = ((rjy — ri)^). (Other 
moments yield similar behavior [3].) We consider the effective bond length b{x,N) = 



Re{x, N) I y/N — 1 to compare the measured chain size with the ideal chain behavior which is 



commonly taken as granted 
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l55l | and which is the basis of our perturbation calculation. 
We use the notation b{x) = lim7v-»oo b{x, N) for the effective bond length of asymptotically 
long chains \^^- The effective bond lengths b{x,N) for = 64 and = 2048 and the 
asymptotic limit b{x) — obtained by extrapolation as described below — are presented in 
Fig. [9] as functions of x. Obviously, b{x, N) l{x = 0) for all N in the small-x limit. 
b{x, N) increases then in the intermediate x- window before it levels off at x ~ 10. The 
swelling due to the excluded volume interaction is the stronger the larger the chain length, 
i.e. b{x, N) increases monotonously with A^. This swelling therefore cannot be attributed 
to a /oca/ persistence length as described, e.g., by the freely-rotating chain model. 

The chain length effects can be seen better in Fig. [TO] where we have plotted b{x, N) for 



several penalties a; as a function of t = 1/y/N — 1. The choice of the horizontal axis is 
motivated by Eq. ([7]) suggesting the linear relation 

b'^ix, N) ^ b'^ix) (1 - c{x)cs{x)t) (27) 

for u = N/g ^ 1 with Cs(x) = A/24/7r^/ pb{x)^ being the swelling coefficient defined in 
Sec. [H c{x) is an additional numerical prefactor of order unity which has been introduced in 



agreement with Eq. (19) of Ref. 12|]. The reason for this coefficient is that the corrections to 
Gaussian behavior differ slightly for internal chain segments [as described by Eq. (I7j)] and the 
;otal chain size which is characterized in Fig.[TOl It has been shown that c — > 1.59 for large 



12 1 . However, since this value corresponds to the limit of a very slowly converging integral 
12l | it is better to use Eq. fl27|) as a two-parameter fit for b{x) and c{x) and to crosscheck then 
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whether the fitted c is of order unity. As shown in the figure for three overlap penalties, this 
method can be used reasonably for overlap penalties as low as a; ~ 0.1, albeit with decreasing 
X systematically underestimating the "true" 6- values indicated in Table [ITl Please note that 
N/g ^ 400 for x = 0.1 and = 8192. Chains with > 8192 would be required to use this 
method for even smaller x. In this limit it is better to use as a simple first step the value 
b{x, N = 8192) of the largest chain length simulated as a (rather reasonable) lower bound 
for b{x). 

Segment size R{s). As we have already stressed in Ref. 12|], it is technically better to 
extrapolate for the effective bond length b using the distribution R{s) of the mean-squared 
size of segments of curvilinear arc-length s defined in Sec. HI It follows from c > 1 that 
the total chain ratio b'^{x, N) converges less rapidly to the asymptotic Gaussian behavior as 
R{sY/s. (See Fig. 4 of Ref. Il2j.) The ratio of i?^(s)/s as a function of s is plotted in the 
inset of Fig. [11] for A^ = 2048 and for several x. As can be seen, it increases systematically 
with segment length s. The swelling levels off for large s, but rather gradually. Therefore it 
would not be appropriate to identify the maximum around s ~ A^ as the asymptotic plateau. 
This again would yield an underestimation of b{x). A more precise method to obtain b{x) 
uses the predicted correction, Eq. ([7]), to the Gaussian limit. We recommend to plot, using 
double-logarithmic coordinates, 1 — R{sY /b'^s as a function of Cs/y/s and to tune b{x) until 
the data for intermediate chain segments with g <^ s N collapses on the bisection line. 
This one-parameter fit yields good estimates down to x = 0.01 where N/g ^ 40. Since the 
corresponding plot is very similar to Fig. 5 of Ref. jlj], it is not reproduced here. We rather 
show in Fig. [TT] a scaling plot motivated by the key relation, Eq. ([5]), which uses our best 
values of g{x) and b{x) for asymptotically long chains (Table [IT|l . The data collapse on the 
theoretical prediction (bold line) is remarkable, especially considering that u covers seven 
orders in magnitude. The "Fixman limit", Eq. ([8]), for m <^ 1 (dashed line) fits the data 
with the smallest overlap penalty x ~ 0.001 confirming the chosen value of b. The limiting 
behavior for u ^ 1 [Eq. ^^^], characterizing an incompressible melt of thermal blobs of 
length g, is indicated by the dash-dotted line. 

Predicting the effective bond length. Up to now, we have used theory to improve the 
fit of b{x), rather than to predict it from the thermodynamic properties and local model 
features such as the bond length l{x). The increase of the effective bond length for weakly 
interacting and asymptotically long polymer melts has been calculated long ago by Edwards 
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[see Eq. (5.55) of Ref. [39|] or Eq. (11) of Ref. [l2|]. Reformulated using our notations and 
substituting the bare excluded volume parameter v{x) by l/gp Q, Q, 47] his result reads 



I 1 + I with = (28) 



vr 



where is the bond length of the unperturbed reference chain of the calculation and Gz{K, g) 
the relevant Ginzburg parameter quantifying the strength of the interaction acting on a chain 
segment of length s = g (see Eq. (6) of Ref. 3|). Since G^ becomes small for large com- 
pressibilities g{x), one expects good agreement with our data for small x. The question is 
now what actually might be the best reference bond length to allow a prediction over the 
broadest possible x-range. The simplest choice to associate bj- with the bond length l{x) 
yields the dash-dotted line indicated in Fig. [91 As can be seen, this choice of b^ allows a 
reasonable prediction only up to x ~ 0.01. The predictive power of Eq. (!28|) can be consid- 
erably improved over nearly two decades up to x ~ 1 if one applies the formula iteratively 
starting with b^ = I and using the effective bond length obtained as input for the Ginzburg 
parameter {b bj) in the next step. This recursion relation converges rapidly as shown by 
the bold line indicated in Fig. [9] obtained after 20 iterations. This iterative renormalization 
of the bond length of the reference chain and the associated Ginzburg parameter has been 



suggested by Muthukumar and Edwards 



67|. Essentially the same result is obtained up to 



X ~ 1 if one sets directly b^ = b using the measured effective bond length (not shown), i.e. 
these values correspond to the fix-point solution of Eq. ( l28l) . The Ginzburg parameters G^ 
obtained using the measured b{x) are listed in Table HTl Note that < 0.34 for x < 1 where 
Eq. ( l28l) fits our data nicely. The fix-point solution of Eq. ( l28i) does not capture correctly 
the leveling off of 6(x) setting in above x ~ 1. Since the Ginzburg parameter becomes there 
of order one (Table HTl). this is to be expected. In summary, we have shown that the iteration 
of Eq. (|28|) allows a good prediction for 6(x) for x < 1 such that G^ <^ 1. If reliable values 
for the compressibility g{x) are available [by means of the extrapolation method implied by 
Eq. ( l23l) ]. this is the method of choice if one cannot afford to simulate very long chains. 



G. Bond-bond correlation function 



Motivation and theoretical prediction. We return now to the deviations from Flory's 
ideality hypothesis predicted in Eq. ([5]) for the mean-square segment size B?{s). As we have 
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seen above (Fig. [TT]), this property requires to substract a large Gaussian contribution b'^s 
from the measured -R^(s) to demonstrate the existence and the scahng of the deviations. 
Unfortunately, this requires as a first step the precise determination of the effective bond 
length b{x) for asymptotically long chains which might not always be available. Indeed we 
have used in the preceeding Sec. IIIIFI the fact that the scaling of Eq. critically depends 
on this accurate value to improve the estimation of the effective bond length b{x) for asymp- 
totically long chains. Hence, it would be nice to demonstrate directly the scaling implied 
by our key prediction without any tuneable parameter. The trick to achieve this is similar 
to our demonstration of the density fluctuation contributions to the free energy, Eq. ([21), 
presented in Sec. IIIIBI We consider the curvature of -R^(s), i.e. its second derivative with 
respect to s, to eliminate the large Gaussian contribution. In principle this can be achieved 
by fitting -R^(s) by a sufficiently high polynomial whose second derivative with respect to 
s then is compared to the theory. Following jo], [l^ we use a more direct numerical route 
where we compute the well-known bond-bond correlation function P{s) = ijrn=n+s ■ In) /^^ 
with li = rj+i — Vi denoting the bond vector between two adjacent monomers i and i + 1 
and the mean-square bond length (Sec. IIII El) . The average is performed as before over 
all chains and all pairs of monomers {n, m + s) possible in a chain of length A^. We use this 
definition rather than the more common first Legendre polynomial (e„ ■ e^) since it allows 
to relate the bond-bond correlation function to the segment size by 

This formula is obtained from (Z„ ■ 1^) ~ {dnTn ■ dmfm) = —dndm {{fn — fmY) /2. Using 
Eq. (|29l) the key prediction, Eq. ([5]), implies for the bond-bond correlation function 



where we have introduced the coefficient cp = Cs(6//)^/8. Eq. (!30l) corresponds to the 
limiting behavior 

ns) « (31) 

for small reduced arc- lengths m <^ 1. The explicit compressibility dependence drops out in 
the opposite limit (m 3> 1) where the bond-bond correlation function becomes 

P{s) ^ c^ls^l\ (32) 
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in agreement with Eq. ([7]). Please note that cp depends imphcitly on the compressibihty. 
(Obviously, both asymptotic behaviors could have been obtained directly from the corre- 
sponding limits for -R^(s), Eqs. ([HD and ([7]).) 

Numerical confirmation. The bond-bond correlation function P{s) for different overlap 
penalties x is presented in Fig. [T^] for chains of length = 2048. As can be seen from 
the unsealed data shown in the inset, P{s) approaches a power law with exponent = 1/2 
(dashed line) in the limit of weak overlap penalties in agreement with Eq. fl3Tl) . For x > 1 our 
data is compatible with an exponent c<j = 3/2 (dash-dotted line) as suggested by Eq. (!32l) . 
Hence, we have demonstrated without any tunable parameter that Flory's ideality hypothesis 
is systematically violated for all segment lengths s and all overlap penalties x. 

We consider now the prefactors and the scaling with x. As suggested by Eq. fl5U]) . the 
main figure presents P{s) / {cp / g{xY^'^) as a function of the reduced arc-length u = s/g{x) 
using the dimensionless compressibilities g{x) and effective bond lengths b{x) from Table HTl 
The data collapse is remarkable as long as 1 ^ s ^ A^. The relation Eq. (l30l) is indicated 
by the bold line; it is in perfect agreement with the simulation data 68|. The asymptotic 
power law behavior with exponents u = 1/2 for m ^ 1 and u = 3/2 for m ^ 1 is shown 
by the dashed and dash-dotted lines, respectively. As predicted by Eq. fl32l) . one recovers 
the power law P{s) = c-p/s^^'^ — already observed for incompressible melts j^, [l^ — for 
scales larger than the thermal blob irrespective of the blob size g. This demonstrates that 
the exponent = 3/2 is not due to local physics on the monomer scale, since ioi s ^ g ^ 1 
distances much larger than the monomer or even the thermal blob are probed. 



IV. CONCLUSION 



Thermodynamic properties of a BFM version with finite overlap penalty. In this paper 
we have discussed a generalization of the standard bond-fluctuation model (BFM) where 
the monomers may overlap subject to a finite energy penalty e (Fig. [1]). This allows us to 
switch on systematically the excluded volume interaction between the monomers as sug- 
gested by perturbation theory [33] and to tune the density fluctuations of the solution at 
constant monomer density. In this study we have focused on dense polymer melts contain- 
ing flexible linear chains which are athermal apart from the finite overlap penalty. The 
central thermodynamic parameter characterizing these systems is the excess part of the 
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dimensionless compressibility g = Tn^p of the solution which has been obtained directly 
from the low-wavevector limit of the static structure factor (Figs. [7] and [8]) . Scanning the 
overlap penalty (or, equivalently, the temperature T [48|) from x = e/T = oo (no overlap) 
down to X = 0.0001 leads to a variation of g{x) over four orders of magnitude (Fig. Ej). 
This allows for a systematic study of the thermodynamic properties (Figs. [2J1SD and the 
intrachain configurational statistics (Figs. [9] fT2l) . Particular attention has been paid to the 
thermodynamic properties of weakly interacting melts (a; ^ 1). We have verified that our 
results are consistent with the free energy, Eq. ([2]), postulated in agreement with Edwards 



391] . The main result of this part of our study is that we have been able to demonstrate the 



density fluctuation contribution to the free energy induced by the chain connectivity from 
the scaling of the specific heat cy with respect to overlap penalty x (Fig. Hj). 

Intrachain conformational properties: Violation of Flory's ideality hypothesis. The 
broad variation of g{x) puts us into a position to test the recently proposed Eq. ([5]) pre- 
dicted by perturbation theory 13|, l53|] describing the systematic swelling of chain segments 



as function of the segment size s and the compressibility g{x). As outlined in the Intro- 
duction, this key relation suggests that the repulsive interactions between chain segments in 
the same chain are not fully screened at variance to Flory's ideality hypothesis for polymers 
in dense melts 40| . The violation of the ideality hypothesis is demonstrated in Fig. [11] for 



the mean-square segment size -R^(s) and in Fig. [12] for the bond-bond-correlation function 
P{s). We show that data obtained for systems with very different compressibilities g{x) 
can be superimposed on the predicted master curves if plotted as a function of the reduced 
arc-length u = s/g. The scaling of R{s) allows a precise determination of an important 
intrachain property, the effective bond length b{g) for asymptotically long chains (Fig. [9]). 
These values cornpare well for x < 1 with the fix-point solution of the recursion relation, 
Eq. ([28]) 12], [39], |67| . The bond-bond correlation function P{s) being the second derivative 
of -R^(s) with respect to s allows an even more direct test of the predicted deviations. The 
reason is that the large Gaussian contribution b'^s, which must be subtracted from -R^(s) 
(see the vertical axis of Fig. [TTl) . drops out due to the differentiation. In contrast to Flory's 
hypothesis, P{s) does not vanish rapidly on scales corresponding to the local persistence 
length (Fig. [T2]) . In perfect agreement with theory 68|, the scaling plot shows two power 
law regimes characterized by exponents u = 1/2 for small (u ^ 1) and uj = 3/2 for large 
(m ^ 1) reduced arc- length. 
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The central result of this study is that even for polymer melts with finite overlap penalty 
excluded volume interactions are not fully screened. If distances smaller than the thermal 
blob size are probed the chains are swollen according to the standard Fixman parameter 
expansion. More importantly, even on distances larger than the thermal blob size {s/g ^ 1) 
deviations from ideal chain behavior are found. Interestingly, in this limit the explicit 
cornpressibility dependence drops out and the relations established for incompressible melts 
[sl, [12I are recovered. This shows that soft melts behave on large scales as incompressible 
packings of blobs. 

Outlook. Since the presented soft BFM is fully ergodic (in contrast to the classical 
BFM) and very efficient, it may be an interesting alternative to vario us pop ular coarse- 



69|, 



grained simulation approaches using soft effective interaction parameters 
The presented model is part of a broader attempt to describe systematica. 

iuctuations in dense polymer systems, both for static 



m 



71 



72 



73|. 



correlated density 
dynamical |l7, 76 



42 



ly the effects of 



43 



74 



75] and 



771] properties. This also involves the comparison with (off-lattice) molec- 



ular dynamics simulation using a standard bead-spring model which is discussed elsewhere 

" n n □ 

131 . 1781 . 179| . An important unresolved question is for instance whether recently pre- 



12 



dieted long-range repulsive forces of van der Waals type ( "Anti-Casimir effect" ) [43|, |75| can 
be demonstrated numerically from specific non-analytic deviations from the RPA formula, 
Eq. fl25l) . at intermediate overlap strengths (x ~ 1). In order to do this, we are currently 
improving the statistics of our data. 

In this paper we have discussed only static properties of the soft BFM. Similar scaling 
behavior also has been obtained for the static Rouse mode correlation function which displays 
systematic deviations from the scaling expected for ideal chains We currently are 

working out how these deviations may infiuence the dynamics for polymer chains without 
topological constraints. (These constraints can be switched off even for x = 00 by using the 
"L26" local moves described in Sec. [Tll) Conceptually important issues can be addressed 
if the (artifical) slithering-snake dynamics is analysed and compared to predictions of the 



76, 



771 ] for real, although 



"activated-reptation dynamics" hypothesis suggested by Semenov 
extremely long polymer chains. If no overlap is allowed, the slithering-snake dynamics is 
mown to show anomalous curvilinear diffusion and correlated motion of neighboring snakes 



171 ]. Since for a; = 00 the lattice might influence the results, it is important to verify 



if qualitative similar behavior is also found for soft BFM melts with thermal blobs much 



26 



larger than the local monomer scale and how the anomalous curvilinear diffusion changes 
with compressibility. 
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TABLE I: Various properties for soft BFM beads {N = 1) at monomer number density p = 0.5/8 
(corresponding to a volume fraction 0.5) and linear box size L = 256 as a function of the reduced 
overlap strength x = e/T. The limit x = oo corresponds to the classical BFM without monomer 
overlap, the limit x = to non-interacting monomers. Indicated are the mean energy per bead e, 
the specific heat cy per bead, the excess part of the chemical potential fiex/T, and the dimensionless 
compressibility g{x,N = 1). Within statistical accuracy we obtain below x 0.1 the ideal gas 
compressibility, 5 ~ 1, and above x ^ 10 the compressibility for a melt without monomer overlap. 
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TABLE II: Various properties for asymptotically long BFM chains at number density p = 0.5/8 
as a function of a; = e/T. Apart from the properties already presented in Table [I] for beads we 
indicate here the intrachain self-energy Cseif, the root-mean-square bond length / = (|Z^^ , the 
effective bond length b, the mean angle {6) and the mean cosine (cos(0)) = (e„ • e^+i) of two 
subsequent bonds with e„ = ln/\ln\ being the normalized bond vector, the swelling coefficient 
Cs = -s/^A/tt^ / pb"^ , and the Ginzburg parameter Gz = 1/ ^b^p. The excess part of the chemical 
potential of a chain is given in units of the chain length N (column 4). The effective bond length 
b{x) has been obtained using an extrapolation scheme implied by Eq. ([7]) and discussed in Sec. lIIIFl 
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TABLE III: Various properties for BFM melts of number density p = 0.5/8 at overlap strength x = 
e/T = 1 as a function of chain length N. For small chains the overlap energy e and its fluctuation cy 
increase slightly while the chemical potential per bead decreases. The compressibility g{x = 1,N) 
becomes chain length independent for N > 64. The chain length dependence visible for small N is 
described by Eq. ([23]) . i.e. the data is consistent with an excess compressibility Qexix) 2.4A^'^ for 
all N. The last column indicates the rescaled end-to-end distance b{x = 1,N) = Rc{N) / {N — 1)^^'^ 
which approaches the effective bond length 6(x) ~ 3.13 of asymptotically long chains monotonously 
from below, just as for classical BFM melts [12]. Interestingly, b{x = 1,N) has not yet reached 
the asymptotic limit b{x) even for N = 8192 albeit all other quantities indicated can be regarded 
(within statistical accuracy) as independent of chain length above ~ 256. 
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FIG. 1: Sketch of the bond- fluctuation model (BFM) with finite excluded volume penalty. The 
BFM algorithm represents monomers by cubes of length a on a simple cubic lattice (of lattice 
constant a) which are connected by a set of allowed bond vecto rs g iven by Eq. Q. Two short 



37l | assumes that all lattice sites 



chains of length = 3 are shown. The classical BFM model |l|, |3|, 
are at most occupied once. We relax this constraint and penalize double occupation by a finite 
interaction energy e which has to be paid if two cubes totally overlap. A corresponding fraction of 
the energy penality is associated with a partial monomer overlap, as sketched in the figure for two 
cube corners occupying the same lattice site. Varying systematically the ratio x = e/T al 



ows us 



to put to a test various theoretical results obtained by perturbation calculation 8 
flexible polymer chains in the melt. 
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FIG. 2: Reduced mean overlap energy per monomer y = e/e as a function of the overlap penalty 
X = e/T for several chain lengths N as indicated. The energy decreases monotonously with 
increasing x. The decay becomes Arrhenius-like for x ^ 10 (dash-dotted line). The dashed line 
indicates the energy predicted from the second virial of soft BFM beads, Eq. (jl5p . The main figure 
demonstrates the weak chain length dependence on logarithmic scales, especially for strong excluded 
volume interactions (x > 1). Inset: Same data plotted with linear vertical axis emphasizing the 
higher mean energy for long polymers {N > 64) for x ^ 1 caused by a self-energy contribution 
Cseif/e ~ 0.18. The self-energies are indicated by the triangles. The bold line shows the temperature 
dependence predicted by Eq. pT|). 
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FIG. 3: Reduced mean energy e/e (spheres) and self-energy Cseif/s (triangles) as functions of the 
number density p for N = 8192, L = 512 and x = 0.001. As shown by the dashed line, e{p) 
is a superposition of the mean field energy 4p and the (essentially) constant self-energy Cseif/e ~ 
O.lSN^x^ p^ . Inset: e/e — Ap as a function of chain length N — \ for our reference density 
p = 0.5/8 and for a single chain {p = 0). The linear slope (bold line) is expected from the return 
probability of Gaussian chains. 
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X=8/T 

FIG. 4: Specific heat per bead cy vs. x for chain length N as indicated. The dashed hne indicates 
the energy fluctuations predicted from the second virial, Eq. (jlSp , which fits nicely the data of soft 
BFM beads {N = 1) over six decades. While the chain length appears not to matter for strong 
excluded volume interactions, the energy fluctuations are found to increase strongly with N for 
X <^ 1. For short chains we observe cy ~ pN^^'^x'^ as can be seen for = 16 (thin solid line). The 
chain length effect saturates for long chains where cy ~ p^^'^x^^'^N^ (bold line) in agreement with 
Eq. ([H]). Inset: cy/Cp^/^x^/^) function of the reduced chain length u = N/g{x) with g(x) 
being the dimensionless compressibility (Table [Tl]) . 
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FIG. 5: Excess chemical potential y = ficx/TN as a function of the inverse temperature x = e/T. 
Increasing linearly (dashed line) for small x it levels off for large x ^ 1 (dash-dotted line). The 
dashed line indicates the simple second virial approximation y ~ v{x)p for unconnected beads, 
fitting successfully the data below x ~ 1. The bold line corresponds to the high temperature 
prediction Eq. (|20p taking into account the density fluctuation contribution induced by the chain 
connectivity. Inset: The chemical potential has been obtained by thermodynamic integration over 
the excluded volume interaction of an inserted ghost chain generalizing the method suggested in 
Ref. [5|. The mean number of lattice sites where monomers and ghost monomers overlap, (Ngg), 
is presented for N = 2048 as a function of A = exp(— Ssg/ST) for x = 3 and x = 100. A power law 
increase of (A'sg) is found for large x (dashed line). 
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FIG. 6: Dimensionless compressibility g{x, N) as a function of x for different chain lengths N 
using the same symbols as in Fig. HI Main panel: Raw data as obtained from the low-wavevector 
limit of the structure factor. Chain length effects become irrelevant for x > 0.1 if > 64 and 
for X > 0.001 if > 2048. The data are compared to the simple second virial approximation 

ity 



1^ 



l/v{x)p (dashed line) which reduces to l/{Sxp) for x <^ 1. As one expects, the compressibi 
levels off for large x and becomes identical to the value g ~ 0.25, known for the classical BFM 
(dash-dotted line). Inset: As suggested by Eq. (p3]) the excess part of the inverse compressibility 
l/g{x,N) — 1/N becomes chain length independent, i.e. the data points for all collapse. The 
master curve indicated by the bold line corresponds to the long chain limit g{x) = limAr-*oo dix, N) 
indicated in Table Ull 
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FIG. 7: Total structure factor G{q) as a function of wavevector q for N = 2048 for different overlap 
penalties x = e/T as indicated. For comparison, we have also included the single chain form 
factor F{q) for x = 0.001. The low- wavevector limit of the structure factor is used to determine 
the dimensionless compressibility g{x,N) [Eq. ([T])]. Only for x < 3 does the structure factor 
decay monotoneously with q as suggested by the RPA formula, Eq. (j25p . G{q) becomes essentially 
constant for smaller temperatures except for wavevectors corresponding to the first sharp diffraction 
peak (called here "Bragg peak"). The box size L = 256 allows only a direct and fair determination 
of g(x,N) for x > 0.1. We have been forced to increase the box size to L = 512 for smaller x as 
may be seen for an example with x = 0.1 (dash-dotted line). As shown by the bold dashed line, 
the RPA formula is used to improve the estimation of g{x, N) for small x. 
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FIG. 8: Rescaled total structure factor G{q)/g{x) as a function of the reduced wavevector Q ^ qC 
for chain length = 8192 and several x < 1 as indicated. The screening length of the thermal 
blob is obtained according to Eq. ^ using the dimensionless compressibility g{x) and the effective 
bond length b{x) from Table [III The bold line compares the data with the approximated RPA, 
Eq. (j26p . If replotted as indicated in the inset the data collapse on the bisection line. Deviations 
from the RPA formula become visible for larger x as shown for x = 1 (crosses). 
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FIG. 9: The (effective) bond length as a function of the reduced overlap penalty x = e/T. The data 
for the root-mean-square bond length l{x) and the effective bond length 6(2;) for asymptotically long 
chains are listed in Tab. [Ul The dash-dotted line indicates the effective bond length as predicted 
by Eq. (j28p assuming = l{x) for the bond length of the reference chain. The bold line shows 
the fix points obtained by iteration of Eq. (j28p using as an input for the Ginzburg parameter the 
effective bond length of the previous iteration step: 6 — > 6r- See the main text for details. 
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FIG. 10: Rescaled mean-square chain end-to-end distance b{x, N)'^ = R'^{N) / (N — 1) as a function 
of t = l/\/N — 1 for different x as indicated. The chains only remain Gaussian on all scales and 
all N for extremely small x. For x > 0.1 one observes b{x,N)'^ to decay linearly in agreement with 
Eq. (j27p . This can be used for a simple two-parameter fit for b{x) as indicated for x = 0.1, 1.0 and 
oo. Note that the coefficient c is slightly above unity as expected from Eq. (19) of Ref. [12l |. 
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FIG. 11: Segment size R{s) for overlap penalty x as indicated for chain length N = 2048. Inset: 
R{s)'^/s as a function of segment length s increases monotonously approaching from below the 
asymptotic limit for large s, i.e. the chains are swollen. Main figure: As suggested by Eq. ([5]), 
the rescaled data (l — i?^(s)/6^(x)s) / [c^{x)/ g{x)^/'^') is plotted as a function of the reduced arc- 
length u = s/g. The data collapse is successful for 1 <^ s <^ N which confirms the values g{x) 
and b{x) for asymtotically long chains (Table [n|) . The bold line shows the full prediction from 
Eq. ([5]). We indicate the limiting behavior for small and large u by the dashed and dash-dotted 
lines, representing respectively Eq. ^ and Eq. ([7]). 
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FIG. 12: Bond-bond correlation function P{s) for different overlap penalties x as indicated in the 
figure. Inset: P{s) as a function of segment length s in log-log coordinates. The data approaches a 
power law behavior, P{s) ~ , with exponent w = 1/2 for small x (dashed line) and uj = 3/2 for 
X > 1 (dash-dotted line). Main panel: Rescaled bond-bond correlation function P{s)/ [cp{g)/g'^^'^'\ 
plotted as a function oi u = s/g as suggested by Eq. ([30|) . For large u, where an incompressible 
packing of thermal blobs is probed, all data collapse onto the dash-dotted line as predicted by 
Eq. (j32p . i.e. P{s) becomes independent of the compressibility g. That this holds not only for the 
classical BFM with x = oo (stars) but also for finite x is the central result of this study. 
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